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ABSTRACT 

We calculate the electronic structure of several atoms and small molecules by direct mini- 
mization of the Self-Interaction Corrected Local Density Approximation (SIC-LDA) func- 
tional. To do this we first derive an expression for the gradient of this functional under 
the constraint that the orbitals be orthogonal and show that previously given expressions 
do not correctly incorporate this constraint. In our atomic calculations the SIC-LDA 
yields total energies, ionization energies and charge densities that are superior to results 
obtained with the Local Density Approximation (LDA). However, for molecules SIC-LDA 
gives bond lengths and reaction energies that are inferior to those obtained from LDA. 
The nonlocal BLYP functional, which we include as a representative GGA functional, 
outperforms both LDA and SIC-LDA for all ground state properties we considered. 



Introduction 



The Local Density Approximation (LDA) |l] has become one of the most popular tools 
for electronic structure calculations. The reason for this is that it gives good accuracy for 
structural properties and is computationally less costly than traditional quantum chem- 
istry methods such as Hartree-Fock, Configuration Interaction and Coupled Cluster meth- 
ods. With the rapid increase in computer power and the development of low complexity 
algorithms, the limits on the system size are being pushed up steadily. However, in many 
cases LDA is not sufficiently accurate. A primary concern is therefore to improve upon the 
accuracy of the LDA approximation at the expense of a moderate increase in computa- 
tional cost. The Generalized Gradient Approximations such as the Becke-Lee- Yang-Parr 
(BLYP) scheme p], H] fall into this category and are now widely used. Numerous other 



schemes to improve upon the LDA scheme can been found in the literature [11|, but very 



few of them have been systematically tested in atomic and molecular calculations. A 
major deficiency of the LDA and also to a lesser extent of the GGA is the fact, that 
there is an unphysical self-interaction in these functionals. To cure this, several years 
ago Perdew and Zunger proposed a scheme |4| where self-interaction terms are subtracted 
out in a straightforward way (SIC-LDA). Even though this scheme is appealing because 
of its conceptional simplicity, it has not been widely used possibly because it leads to 
a numerically more complicated scheme since the potential becomes orbital dependent. 
The minimization of this functional therefore cannot anymore be considered as an eigen- 
value problem in a self-consistent potential and the total energy is not any more invariant 
under unitary transformations among the occupied orbitals. In this paper, we first derive 
the gradient of the SIC-LDA functional, which is necessary for minimization algorithms 
|| [2(| , and then present results that we obtained for atomic and molecular systems. 



Studies of periodic systems using SIC-LDA have mainly concentrated on systems where 
the LDA approximation gives qualitatively wrong results (such as transition metal oxides 
0) and where the electronic structure undergoes qualitative changes in response to chang- 
ing external conditions |7], ||. Most of these calculations involving heavy atoms were done 
with the LMTO method. The gap in insulators is also improved compared with the 
LDA case 0. In the case of atoms and molecules most papers [[U], [12], [L3| investigated 
the properties of excited states of atoms and molecules at the experimental geometry and 
found that the SIC-LDA improves upon LDA. Instead, we concentrate on the ground state 
properties such as the equilibrium geometries of small molecules and find that SIC-LDA 
is less accurate than LDA. 
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The SIC-LDA functional and its gradient 



In all the following formulae we consider the orbitals to be real and the subscripts of 
the orbitals run over all the occupied orbitals. The SIC-LDA functional is given by 



+ i f r p(^m drdrl _ i E f f PfP^) drdr ' 
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where 

Pi(r) = ^(r)^(r), p{r) = Y J Pi{r) 

i 

and / e(p)pdr is the local approximation of the exchange correlation functional. 

We want to obtain its gradient under the constraint that the orbitals be orthonormal. 
Following the ideas outlined by Arias et al. [16] , we consider a more general functional 
that is also defined with respect to nonorthogonal functions First we construct a set of 
orthonormal orbitals by a symmetric Lowdin orthogonalization of the non-orthogonal 
set 

^ = E^ 1/2 ^- (2) 

j 

where Sij =< 9i\9j >= J 9i(r)9j(r)dr is the overlap matrix among the occupied or- 
bitals. The functional we are interested in is just the SIC-LDA functional evaluated for 
the orthonormal orbitals In our actual calculation, for reasons of numerical stability, 
we use orthogonal orbitals. Hence, in our derivations it is necessary to consider only in- 
finitesimally nonorthogonal orbitals . Then S~ 1/2 — (I + (S — I))~ 1/2 « / - (1/2) (5 - I) 
and Eq. (2) simplifies to 

*i = E(f «w - (3) 

The gradient of the total functional is then obtained by applying the chain rule: 

-gg- = W dE d -Ml dr > (a) 

where each part of Eq. (4) can easily be calculated. Denoting the unconstrained gradient 
by dj(r) we obtain: 

1 f)W /I \ 
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where the orbital dependent Hamiltonian Hj is 

1 
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The second part of Eq. (4) gives 



dtfi(r) 



^(r-r')-^(r-r')-^E * 
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E(^(r) + ^(r))^(r') 



(6) 



In the last transformation step, we have used the fact that we calculate the derivative 
for a set of orthonormal orbitals and therefore S = I. In order to take account of the 
orthogonality constraint, we are of course allowed to put S = I only after calculating the 
derivative. Finally we obtain the gradient expression 



1 dE 



di(r) 



2dUr) — " ^(/^(rOf)*i(r) 
- ^E(/^( r/ )^(r')rfr')^(r). 
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The above gradient expression is different from what is found in the literature [12, |15 , [7] . 
1 dE 



2d* At) 



^W-E(/^( r, )^(r')^') *,-(r) 



In numerical applications, it was apparent that the gradient expression in Eq. |8] does not 
lead to the correct minimum []15] of the functional. However, nobody seems to have drawn 
the logical conclusion that Eq. |8] is not the correct gradient of the LDA-SIC functional. 
Instead people added a second minimization step which is based on a relation derived by 
Pederson et al HlO 

< VilHM >=< H^l^j > . (9) 
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This relation follows immediately by considering infinitesimal unitary transformations 
among the occupied orbitals and putting the gradient of the SIC-LDA functional with 
respect to these transformations equal to zero. We see that Eqs. || and || together are 
equivalent to Eq. 0. These difficulties appear only in the SIC-LDA case; in the LDA 
case the potential is not orbital dependent and therefore the two gradient expressions are 
identical. 

The results of the following two sections were obtained by direct minimization of the 
SIC-LDA functional, using the correct gradient [F[ The DIIS method []B5| was used as 
a convergence accelerator. In the molecular case the gradient was preconditioned using 
the scheme by Teter et al. [20|, in the atomic case the operator Im{^~) was used as 
a preconditioner ]2T| , where z is a suitably chosen complex energy. All the calculations 
were done for a spin unpolarized system where the spin up orbital is required to have the 
same spatial form as the matching spin down electron orbital, \l/2i-i(r) = ^2i(i*) 

To calculate the equilibrium geometries of small molecules, we relaxed the atoms in 
the direction of the forces until the forces vanished. The forces in the SIC-LDA scheme 
are given by the Hellman-Feynman theorem. This might not be quite obvious, since in the 
usual derivation of the Hellman-Feynman theorem one takes advantage of the fact, that 
the orbitals are eigenf unctions of the self-consistent Hamiltonian, which is not the case 
in the SIC-LDA scheme. We therefore give in the following a derivation of the Hellman- 
Feynman theorem which does not require the orbitals to be eigenfunctions, but uses only 
the fact that the orbitals minimize the total energy for some fixed positions of the nuclei. 

Let us consider a set of orbitals \&j(r, R') which are the solutions of the SIC-LDA 
equations for a molecule, whose atomic positions are given by R'. We consider R' to be 
a super vector, i.e. for a system comprising N atoms it will have 3N components. By 
construction these orbitals are orthonormal for all values of R. Let us now consider the 
SIC-LDA total energy E tot [^i(r, R'); R] for a set of atomic positions R. The dependence 
on R stems from the fact that the external potential V ext depends on the atomic positions 
R. Obviously the functional will be minimized if ^i(r, R') = \&j(r,R) and its gradient 
with respect to R' will vanish at that point. 



<9£ tot [^(r,R);R] 



<9R 



R'=R 



(10) 



Let us now calculate the force acting on the nuclei which is given by 



<9£ tot [^(r, R); R] <9£ tot [^(r, R); R" 



OR 



OR" 



+ 



dE tot [%{r,R');R] 



R"=R 



OR' 



(11) 



R'=R 



The first term on the right side of Eq. (11) takes into account that the external potential 
V ext depends on the atomic positions R but freezes the dependence of the orbitals on R. 
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It is given by 

W^( r ,R)-^^(r,R)dr (12) 

The second term freezes the dependence on the external potential but takes into account 
the dependence of the orbitals . This second term is however equal to the left hand side 
of Eq. (10) and therefore vanishes. So we are just left with the first term which can be 
immediately transformed to the usual Hellman-Feynman form. 



Atomic results 

In their original paper, Perdew and Zunger [|J do several atomic calculations. They 
however did not minimize the SIC-LDA functional under the orthogonality constraint, 
but instead employed the eigen orbitals of the orbital dependent potential. Since the 
orbitals are the eigenfunctions for different Hamiltonians, they are not orthogonal but 
they assume that this will not change their results appreciably. We therefore repeated 
some of the atomic calculations to check this assumption and in fact find it to be true. 

The orbitals that one obtains from a minimization solution differ in two important 
aspects from the atomic orbitals obtained in a non-orbital dependent potential. First 
their nodal structure is different as shown in Fig.l. All orbitals (even the Is core states) 
have the same number of nodes. Second all the orbitals decay with the same exponent (see 
Fig. 2). A set of orbitals with a behaviour appropriate for a local potential can however 
be obtained by forming the new linear combinations 

= ( 13 ) 

j 

where C is the matrix containing the eigenvectors of the (hermitian) matrix < ^>i\H^ j >. 

In Table [I], we compare our total energies for several closed shell atoms with the ones 
from Perdew and Zunger. In spite of the fact that the minimizing orbitals have a 
quite different behaviour from the approximate orbitals calculated by Perdew and Zunger 
(which show the behaviour of eigenorbitals in a local potential and look therefore similar 
to the orbitals $j of Fig.l), the total energies are nevertheless very similar. In fact the 
small differences in the total energy in the columns labelled SIC-LDA and PZ are not 
due to their approximate solution method but are instead probably due to the use of an 
insufficiently dense grid in Ref. [Q, since we find the same differences in comparing the 
LDA results too. 

In Table [TJ, we also give the results obtained from the Perdew-Zunger (PZ-LDA) and 
the Perdew- Wang '92 (PW92-LDA) parametrizations, the Becke-Lee- Yang-Parr GGA 
(BLYP) and Hartree Fock (HF). Of the two LDA functionals, PZ-LDA and PW92-LDA, 
the latter is the more accurate parametrization of the correlation energy of a homogeneous 
electron gas and yields slightly more accurate total energies. We include the former only 
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Comparison of Minimizing and Eigen Orbitals of Ar 




-2.0 1 ' ' ' ' ' 1 

0.0 0.5 1.0 1.5 2.0 2.5 3.0 

r (a ) 



Figure 1: Is, 2s, and 3s orbitals of Argon. Solid line: The orbitals ^ that minimize the 
SIC-LDA functional, dashed line: the linear combinations $ of Eq.13. 




Figure 2: Same as Fig.l but on a logarithmic scale. The 3s orbitals \I/ and $ are nearly 
indistinguishable on this plot beyond r — 1. 
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to facilitate comparison with Ref. [|]]. As observed in earlier papers LDA yields total 
energies that are not sufficiently negative while SIC-LDA gives too deep total energies. 
The absolute value of the error is significantly smaller in SIC-LDA than in LDA, but not 
as small as for BLYP. 

In Table [| we compare the highest occupied eigenvalues of closed shell atoms. Here 
the SIC-LDA functional outperforms all the others. 

In Figs. 3, 4 and 5 we show the error in the self consistent charge densities from 
the different methods. The SIC-LDA density is somewhat more accurate than the LDA 
density except in a region around 0.3 ao for Neon. 



Table 1: Total energies of closed shell atoms in eV. Our results using direct minimization 
and the PZ parametrization (SIC-LDA) are compared with the results obtained by Perdew 
and Zunger (PZ) for the same functional using an approximate solution method. The next 
two columns show the LDA results with the Perdew-Zunger (PZ-LDA) and the Perdew- 
Wang '92 (PW92-LDA) parametrizations. The final three columns are the Becke-Lee- 
Yang-Parr GGA (BLYP), Hartree Fock (HF) and the exact value. The exact values are 
from Ref. [17]. We have used a factor of 27.2112 eV/Hartree for converting most of the 
energies but for consistency with Ref. [|J a factor of 27.21 eV/Hartree for converting the 
SIC-LDA energies. 





SIC-LDA 


PZ 


PZ-LDA 


PW92-LDA 


BLYP 


HF 


exact 


He 


-79.4 


-79.4 


-77.1 


-77.1 


-79.1 


-77.9 


-79.0 


Be 


-399.9 


-399.8 


-393.1 


-393.1 


-398.9 


-396.6 


-399.1 


Ne 


-3517.9 


-3517.6 


-3489.1 


-3489.3 


-3509.2 


-3497.9 


-3508.6 


Ar 


-14378.8 


-14378.3 


-14307.9 


-14311.6 


-14352.4 


-14335.4 


-14354.6 



Table 2: Highest occupied eigenvalues of closed shell atoms in eV. Our results using di- 
rect minimization and the PZ parametrization (SIC-LDA) are compared with the results 
obtained by Perdew and Zunger (PZ) for the same functional using an approximate so- 
lution method. The next two columns show the LDA results with the Perdew-Zunger 
(PZ-LDA) and the Perdew- Wang '92 (PW92-LDA) parametrizations. The final three 
columns are the Becke-Lee- Yang-Parr GGA (BLYP), Hartree Fock (HF) and experiment. 
The experimental values are from Ref. ||18|| . 





SIC-LDA 


PZ 


PZ-LDA 


PW92-LDA 


BLYP 


HF 


exper 


He 


-25.8 


-25.8 


-15.5 


-15.5 


-15.8 


-25.0 


-24.6 


Be 


-9.1 




-5.6 


-5.6 


-5.4 


-8.4 


-9.3 


Ne 


-22.8 


-22.9 


-13.5 


-13.5 


-13.2 


-23.1 


-21.6 


Ar 


-15.9 


-15.8 


-10.4 


-10.4 


-10.0 


-16.1 


-15.8 
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Helium Charge Density Differences 




r (a ) 



Figure 3: Comparison of the charge density of He obtained by different methods with 
the quasi exact charge density from a Hylleras-type calculation. 
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Figure 4: Comparison of the charge density of Be obtained by different methods with 
the quasi exact charge density from a Quantum Monte Carlo calculation. 
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r (a ) 



Figure 5: Comparison of the charge density of Ne obtained by different methods with 
the quasi exact charge density from a Quantum Monte Carlo calculation. 



Molecular results 

Using pseudo-potentials and plane waves as basis set, we did a series of molecular calcu- 
lations. We examined the bond lengths of several small molecules and the energy released 
in a chemical reaction. We generated SIC-LDA pseudo-potentials according to the pro- 



cedure described in reference 2^ where however the eigenvalues and charge distributions 



of the reference configuration were the ones from an atomic SIC-LDA calculation. As 



is to be expected from the concept of pseudo-potentials representing physical ions JL4 
these pseudo-potentials are very similar to the pure LDA pseudo-potentials and substi- 
tuting a LDA pseudo-potential in a SIC-LDA calculation has only a very minor effect. 
For consistency, we used however always the SIC-LDA pseudo-potential. With this kind 
of pseudo-potential it is possible to get the correct LDA bond length to within a few 



thousands of a Bohr for first row molecules WM. We expect the same accuracy for the 



SIC-LDA pseudo-potentials. A very attractive feature of the SIC-LDA scheme is that the 
minimizing orbitals can usually easily be interpreted in physical terms. They represent 
either bonds or lone electron pairs. In the case of the CH4 molecule for instance, one 
obtains 4 localized orbitals (each containing a spin up and spin down electron), which are 
centered on the 4 lines linking the 4 hydrogens with the central carbon and representing 
therefore bonds. In the case of the H2O molecule one obtains again 4 localized orbitals 
in nearly tetragonal positions. In this case however only the two sitting on the two lines 
linking the O with the two H represent bonds, the other two which are in the half-space 
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not containing any H are lone electron pairs. In the case of double or triple bonds, the 
minimizing functions are banana shaped localized functions distributed around the line 
linking the two atoms. We found that in the case of a Si crystal also, the minimizing or- 
bitals are bond-centered Wannier type functions. In some cases such as the CO molecule 
we found two very close minima, the lowest one corresponding to a triple bond, the other 
one to a single bond. The charge density was however very similar in both cases. Unfor- 
tunately, as can be seen from table [J[ the bond lengths we obtain for small molecules are 
systematically too short, and the error was significantly larger than the error of the LDA 
bond length. The error for triple bonds are particularly large. 



Table 3: Comparison of the LDA and SIC-LDA bond lengths (a.u.) for several small 
molecules. The experimental bond lengths (exper.) and the differences between the 
theoretical and experimental bond lengths are given. 





exper. 


SIC error 


LDA error 


BLYP error 


H 2 


1.401 


-.03 


.05 


.01 


CH 4 


2.052 


-.05 


.02 


.02 


C 2 H 2 (CH) 


2.005 


-.05 


.03 


-.01 


C 2 H 2 (CC) 


2.274 


-.09 


-.01 


.00 


NH 3 


1.912 


-.05 


.02 


.02 


H 2 


1.809 


-.05 


.02 


.03 


BH 


2.329 


-.09 


.04 


.01 


LiH 


3.015 


-.06 


.01 


.00 


N 2 


2.074 


-.09 


-.01 


.01 


CO 


2.132 


-.10 


-.00 


.02 


co 2 


2.192 


-.09 


.00 


.02 


RMS deviation 




.072 


.024 


.016 



Another interesting quantity in this context is the atomization energy, which is the 
difference between the molecular total energy and the sum of the total energies of its con- 
stituent atoms. To calculate the energy of an open-shell atom one has to decide whether 
one wants to spherically symmetrize the atom by introducing fractional occupation num- 
bers. In the case of GGA is has been found that the non-spherical atom gives better 
results. The theoretical foundation of this empirical observation is however not quite 
clear and it can not be predicted which scheme would give better results in the case of 
SIC-LDA. To check the accuracy of the SIC-LDA total energy differences and avoid at 
the same time the above mentioned symmetry problems we looked therefore at the energy 
released in the chemical reaction 3 H 2 + N 2 — > 2 NH3. All the total energies were cal- 
culated after a full relaxation of the atomic positions within each scheme. In the case of 
SIC-LDA, we also calculated the energy difference using the more accurate LDA geome- 
tries of the molecules. As can be seen from Table 0], the SIC-LDA scheme did worse than 
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the other schemes examined. By far the best results were obtained with the BLYP [[J 
scheme. To check the accuracy of the pseudo-potential plane wave method we calculated 
the LDA result both with this method and with the Gaussian 94 (G94) program 
package, getting extremely close agreement. The BLYP calculation was done with the 
Gaussian 94 software. In all the Gaussian 94 calculations mentioned in this paper we used 
a 6-311G++(3df,3pd) basis set. 



Table 4: The experimental value and theoretical predictions for the energy (eV) released 

in the model chemical reaction 3 H 2 + N 2 — > 2 NH 3 . 

The values are corrected for the zero point energy given in reference [26fl. 



exp. 


LDA (PSP) 


LDA (G94) 


SIC-LDA 


SIC-LDA(LDA geom.) 


BLYP 


.76 


2.1 


2.1 


2.6 


2.9 


.65 



Conclusions 



The SIC-LDA scheme as proposed by Perdew and Zunger does not give ground state 
energies and molecular geometries with sufficient precision. Whereas the atomic results 
are superior to the LDA results, the molecular results are clearly worse. A good GGA 
scheme such as the BLYP scheme outperforms both of them in all test cases considered. 
It is surprising how well the simple LDA scheme works compared to more sophisticated 
schemes that one would expect to give superior results. Charlesworth recently came to 
similar conclusions when he systematically examined several weighted density functional 
schemes |f27|l . All these schemes satisfy sum rules that are generally believed to be 
responsible for the accuracy of the LDA scheme. Thus it might well be that we actually 
do not yet fully understand the true reasons for the success of the LDA scheme. 
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